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We present a self-consistent numerical approach to solve the Gutzwiller variational problem for 
general multi-band models with arbitrary on-site interaction. The proposed method generalizes and 
improves the procedure derived by Deng et a/., Phys. Rev. B. 79 075114 (2009), overcoming the 
restriction to density-density interaction without increasing the complexity of the computational 
algorithm. Our approach drastically reduces the problem of the high-dimensional Gutzwiller mini- 
mization by mapping it to a minimization only in the variational density matrix, in the spirit of the 
Levy and Lieb formulation of DFT. For fixed density the Gutzwiller renormalization matrix is deter- 
mined as a fixpoint of a proper functional, whose evaluation only requires ground-state calculations 
of matrices defined in the Gutzwiller variational space. Furthermore, the proposed method is able 
to account for the symmetries of the variational function in a controlled way, reducing the number 
of variational parameters. After a detailed description of the method we present calculations for 
multi-band Hubbard models with full (rotationally invariant) Hund's rule on-site interaction. Our 
analysis shows that the numerical algorithm is very efficient, stable and easy to implement. For 
these reasons this method is particularly suitable for first principle studies - e.g., in combination 
with DFT - of many complex real materials, where the full intra-atomic interaction is important to 
obtain correct results. 
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I. INTRODUCTION 



In the 60th's Martin Gutzwiller published a series of 
papers 1 3 where he introduced a variational method for 
studying ferromagnetism in transition metals. His bril- 
liant idea was to variationally determine a projected 
wavefunction represented as 



i*g> = n p * i*°> > 



(i) 



R 



where the local operators Vyi improve the non-interacting 
wavefunction |^o) in accordance with the on-site inter- 
action by modifying the weight of local electronic config- 
urations. 

In spite of its simplicity, the average values of any oper- 
ator on \^g) can on ly be computed numerically for realis- 
tic lattice models, e.g., using variational Monte Carlo A£ 
For this reason, Gutzwiller introduced an approximate 
scheme, known as the Gutzwiller approximation, to com- 
pute these average values analytically. Successively, the 
development of dynamical mean field theory (DMFT) 6 
has brought additional insights into the physical meaning 
of the Gutzwiller approximation. In fact, Metzner and 
Vollhardt showed that this approximation is exact in the 
limit of infinite coordination lattices^"— where the single- 
particle self-energy becomes purely local in spaced 

Since its introduction the Gutzwiller wavefunction and 
approximation have proven to be very important tools 
to study strongly correlated systems. The understand- 
ing of many basic concepts, such as the Brinkman-Rice 
scenario for the Mott transition, 11 came originally from 
calculations based on the Gutzwiller method. From the 



computational point of view the list of interesting results 
that have been obtained by means of the Gutzwiller ap- 
proximationi — and its respective generalizations^^ - — 
is impressively long, hence impossible to cite in an ex- 
haustive way. Furthermore, the Gutzwiller approxima- 
tion can be naturally combined with density functional 
theory (DFT) / 6 i 17 applying, e.g., the local density ap- 
proximation (LDA)i 8 for the exchange and correlation. 
LDA+Gutzwiller (LDA+G)^^ has proven to be a pow- 
erful scheme for the study of real strongly-correlated 
metallic materials^ giving a more accurate descrip- 
tion than LDA+U, 21 comparable to LDA+DMFT— for 
ground state properties. Finally, the range of applica- 
tion of the Gutzwiller method has recently been extended 
to out-of-equilibrium calculations, such as electron trans- 
port across quantum dot systems^ 3 - and quench dynamics 
in correlated electron systems*^ 

The Gutzwiller variational method requires a number 
of preliminary technical steps in order to make it re- 
ally flexible and able to cope with systems of interest, 
as the number of variational parameters scales exponen- 
tially with the number of correlated orbitals involved in 
the calculation. This scaling is already problematic for 
transition metal systems with correlated d orbitals if the 
required minimization is performed in a naive way. Based 
on the formalism introduced by Biinemann^ Deng et 
alj^ recently derived a self-consistent numerical method 
that allows to efficiently perform calculations even for 
d-orbital systems. The only limitation of this approach 
is the restriction to density-density type of interactions, 
which is actually due to the employed formalism and does 
not stem from the numerical approach in itself. 

However, in order to properly describe the physics 
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of several strongly correlated materials, the full intra- 
atomic interaction is needed; not only its density-density 
component. The rotational invariant on-site Coulomb 
and exchange interaction is generally modeled in terms of 
the so called Kanamori parameters 2 ^ commonly referred 
to with the symbols U and J. The strength of the two- 
electron spin-exchange interaction is determined by the 
parameter J. In transition metal oxides with partly filled 
d-shells the off-diagonal interactions - exchange coupling, 
spin- flip and pair hopping - are crucial. For example, 
the metallic property of SrV03 can not be reproduced 
from theory without accounting for spin-exchange in- 
teraction^ Several theoretical model studies points in 
the same direction. To mention a few, the transition 
from the paramagnetic to the ferromagnetic phase for 
multiband systems requires finite J ^ the spin-freezing 
transition predicted for multiband systems, which is ex- 
pected to influence the Mott transition^ takes place 
when < J/U < 1/3 and is absent for J = 0, 27 Pre- 
vious Gutzwiller model studies^ 3 - and the present work 
show that the critical U required for the Mott insulator 
transition is substantially reduced when increasing the 
ratio J/U from zero. 

Motivated by the above examples we believe that 
an implementation of the Gutzwiller variational method 
valid for general on-site interactions and, at the same 
time, numerically efficient would constitute an important 
progress. Formal advancements pointing in this direc- 
tion have lately been conceived by several authors j 23 i 28 ~ — 
although these progresses have been hampered by the 
lack of efficient numerical algorithms applicable to the 
most general case. In particular, Fabrizio and collabo- 
rator a 3Q i 32 ~ — derived a mathematical formulation of the 
problem whose complexity is unaffected by the form of 
the on-site interactions and furthermore allows to eas- 
ily incorporate symmetries into the variational function 
from the onset. 

The main goal of this work consists of merging to- 
gether the general formalism developed by Fabrizio and 
collaborators 3 ^ 2 - - — mentioned above and the numerical 
procedure derived by Deng et al^ overcoming the re- 
striction to density-density interaction without increas- 
ing the complexity of the computational algorithm. Fur- 
thermore, the fully rotational invariant on-site interac- 
tion enables us to construct the variational wavefunction 
using the orbital rotational symmetry, which is instead 
broken if the off-diagonal terms are neglected. It will be 
shown that this extra symmetry can be used to reduce 
the number of variational parameters. 

The outline of this paper is as follows. In Sec. [IT] the 
Gutzwiller problem for a general tight binding Hamilto- 
nian is introduced. In Sec. |TTH the employed formula- 
tion of the Gutzwiller metho d 3Q i 32-34 i 36 is summarized. 
In particular, in Sec. IIII CI it is shown that this formu- 
lation provides a natural extension of the formalism of 
Ref. [20, having the same mathematical structure in the 
special case of density-density on-site interactions. In 
Sec. IIVI we discuss the implementation of symmetries of 



the wavefunctions. In Sec. El the numerical procedure 
to minimize the Gutzwiller energy is described in detail. 
In Sec. EH we briefly discuss how the proposed method 
can be adapted to a of LDA+G type of calculation. In 
Sec. EII] we prove the reliability of the method presenting 
a comparison with other methods for the cases of two and 
five orbit als. Furthermore we discuss several technical 
details of the numerical procedure, such as convergence 
properties and computational speed. Finally, Sec. IVIIII 
is devoted to the conclusions. 



II. THE GUTZWILLER METHOD 



Let us consider the general tight binding Hamiltonian 



n 



+ a(3 c ] i 



££[/(R)rr<|r\R)<r',R| 



(2) 



R IT' 

T + H\ oc , 



where cj^ a creates an electron in state a (where a labels 
both the spin a and the orbital a at site R) and |T, R) are 
many-body Fock states expressed in the CR a -basis. These 
states are defined by the occupation numbers n a (T, R) G 
{0, 1} where a runs over integer numbers from 1 to M, 
M being the number of on-site single particle states, 



|r,R) = (4 1 )" l(r,R !. (c RM ) 



uw(r,R) 



10} • (3) 



Thus the number of Fock states is 2 M . The Hermitian 
matrix U (R) represents the local terms, interaction and 
crystal fields, in the cj^ a -basis, i.e., the same basis in 
which T was defined in Eq. ©. This basis will henceforth 
be denoted as the original basis. 

The structure of the Gutzwiller variational function is 
given by Eq. (pQ), where |^o) is an uncorrelated varia- 
tional wavefunction, that satisfies Wick's theorem, and 
Vr is a general operator acting on the local configura- 
tions at site R 



P R = ^A(R)rr'|r\R)(r',R|, 



(4) 



rr / 



where the 2 M x 2 M matrix A(R), assumed to be real in 
this work, contains all the variational parameters needed 
to define the operator Vn- 

In general, average values of operators with respect 
to \S&g) must be computed numerically unless the lat- 
tice has infinite coordination number, in which case they 
can be evaluated analytically if the following equations - 
commonly named Gutzwiller constraints - are satisfied: 



<*o|^P R |*o> 

(W r P r Cr|*o> 



1 

(*o|Cr|* >, 



(5) 
(6) 
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where Cr is the local single-particle density-matrix oper- 
ator with elements cj^Cpo. 

The variational problem to solve amounts to variation- 
ally determine both |^o) and V-r by minimizing the av- 
erage value of the Hamiltonian [Eq.fl2J)] 

£™[P,*o] = (*o\V*HV\* } (7) 
fulfilling Eqs. dSE]), where we have introduced 

v = n p r . (8) 

R 

For a general tight binding model [Eq. (j2j)] this prob- 
lem is complicated for two reasons: (i) |^o) and V are 
not independent variables because of the Gutzwiller con- 
straints [Eqs. (06])], (ii) the number of variational param- 
eters scales exponentially with the number of orbitals. 

III. REFORMULATION OF THE GUTZWILLER 
PROBLEM 

In this section we briefly summarize the reformula- 
tion of the Gutzwiller problem derived in Refs. [33|, [34], 
and we show its formal analogy with the formulation 
of Biinemann and Weber— in the special case of pure 
density-density local interaction. 

A. The mixed-basis representation 

Let us introduce the so-called natural-basis 3 - operators 
d Ra , i.e., the operators such that 

(*ol4a d Rfll*o> = ^ng la =n° i8 (R) Va,/3 (9) 

where < n RQ! < 1 are the eigenvalues of the local den- 
sity matrix 

<*o|cLcr0I*o>=/&0(R). (10) 

Notice that the natural-basis operators are always well 
defined as p°(R) is Hermitian, implying that there always 
exists a unitary transformation U n such that 

4. = E^ c V ( n ) 

Instead of expressing the Gutzwiller projector in terms 
of the original basis as in Eq. (|4]) we adopt the following 
mixed original-natural^ basis form 

P R = ^A(R)m|r,R)(n,R| (12) 

Tn 

where, by assumption, |T, R) are Fock states in the origi- 
nal CR a -basis , while |n, R) are Fock states in the natural 
basis, namely in terms of the -operators. In other 



words a generic state |n, R) is identified by the occupa- 
tion numbers n#(n,R) G {0,1} - with f3 G {1,..,M} - 
and has the explicit expression 

/ . \ni(n,R) / , \n M (n,R) 

KR) = (4i) ••• (4m) |o>. (13) 

For later convenience we adopt the convention that the 
order of the |T, R) and the |n, R) states is the same. For 
instance, if the second T- vector in Eq. (fT2]) is c^cj^ |0) 
then the second n- vector is d\^d\^ |0). 

B. The (j)- matrix 

Let us introduce the uncorrelated occupation- 
probability matrix P°(R)22 with elements 

[P°(R)W = (*olK,R)(n,R||* ) 

= * nn -P n °(R), (14) 

where 

M 

P„°(R) = J] «^) n ' 3(n ' R) ( 1 - n° R/3 ) W "' R) . (15) 

13=1 

We remind that n R/3 are the elements of the diagonal 
density matrix of Eq. (j9]), and denote the occupation 
numbers of the natural states /3. We also introduce the 
matrix representation of the operators and c R ^ 

d R p^ (d R0 ) nni = (n,R\d R0 \n',R) (16) 

c R/3 ^ (c R0 ) rr , = (r,R|c R/3 |r',R>. (17) 

Notice that, if we respect the convention that the order 
of the |T,R) and the |n,R) states is the same, we have 
that 

We now define the matrices A(R) and U(H) with ele- 
ments A Fn (R) [Eq. CL2])] and /7 r r'( R ) [Eq. ©]• Notice 
that A(R) is defined in the original-natural and U(H) 
is defined in the original-original basis. With the above 
definitions, the expectation value of any local observable 
can be calculated as 

(*o| P f O(R)P |*o} = Tr (P°(R)A t (R) O(R) A(R)) , 

(19) 

where 

0^(^ = ^1010, (20) 

and the Gutzwiller constraints [Eqs. (j5]|6])] can be written 

as 

Tr(P°(R)A t (R)A(R))=l (21) 
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Tr(p°(R)At(R)A(R)/t/ /3 )=(* |4 a d R/3 |* >. (22) 

The formalism is further simplified by defining the ma- 
trix 



0(R) = A(R) \/P°(R), (23) 

that was introduced in Ref. [33|. The expectation value of 
any local observable is given by 

(ttol V^6(R)V |* > = Tr (>(R)t O(R) 0(R)) . (24) 
The Gutzwiller constraints take the form 

Tr(0t(R)0(R)) =1, (25) 
Tr (V(R)<KR) /l/^) = <*o| 4^r/3 l*o> 

= ^4, (26) 

and the variational energy [Eq. (|7J)] is, in the Gutzwiller 
approximation, given by 3 - 3 - 

RR' 75 

+ ^Tr(<A(R)t[/(R)0(R)) , (27) 

R 

where 

= EC'W(RV (28) 

K (R) oS . ,20) 
Jn°(R)(l - n°(R)) 

In conclusion, within the formalism summarized in this 
section, the variational energy is a functional of 0(R) 
and l^o)? to be minimized fulfilling the Gutzwiller con- 
straints [Eqs. (I25ll26|) ]. 



C. Diagonal projector as a particular case 

Let us assume that the coefficients of the matrix A that 
define the projector V in Eq. (|4]) are diagonal 

A(R)rr' = $rr> A(R) rr , (30) 

and that the original basis coincides with the natural 
basis 

<^0|4a C R/^0) = ^/3(^0|4c C Rcl^0) 

= <*o|dL^I*o), (31) 
From Eqs. (|30ti3ip we have that 

0(R)y = % 0(R)« , (32) 



and, consequently, the following equation hold for all lo- 
cal operators O(R) 

Tr (0 f (R) <£(R) O(R)) = Tr (0 f (R) O(R) 0(R)) . (33) 

From Eq. ([33]) follows that the equations that charac- 
terize the Gutzwiller problem [Eqs. (|25ti29)) ] can be eval- 
uated in terms of the variational parameters ^ym^(R) 
defined by Biinemann in Ref. [23| 

m)rr = ^/(*o|^t|r,R>(r,R|P|* > 

= Vm r (R). (34) 

A crucial observation in this work is that in the gen- 
eral case considered here, in which Eqs. (|30H3ip are not 
assumed, Eqs. ()25ti29|) are expressed in terms of quadratic 
forms of the matrix elements of 0(R) instead of ^rar(R), 
but in a formally identical way. 

We conclude this section observing that the physical 
density matrix of the system 

p Q/3 (R) = (*o|^ f cLc R ^l*o) 

= Tv(0t(R)/t //3 ^( R) ) (35) 

is not equal to the so called variational density matrix 
n%(K) = (* |dL^I*o> 

= Tr(f(R) </>(*.)&(,) (36) 

when Eqs. (|30H33p do not hold. In the general case the 
distinction between variational density matrix and phys- 
ical density matrix needs to be taken into account. 

IV. SYMMETRIES OF THE VARIATIONAL 
FUNCTION AND MATRIX 

In this section we discuss in detail how to build sym- 
metries in the Gutzwiller variational function. The site 
label R is dropped for simplicity. The procedure dis- 
cussed here extends the method discussed in Ref. [34] to 
general point symmetry groups. The problem amounts 
to define the form of the ^-matrix such that \^g) is in- 
variant under the action of a matrix representation of a 
group G in the many-body R-local space. The transfor- 
mation law 

0<£$ _1 = I>0a(0)4 VjeG 

ff|0) = |0) (37) 

defines a representation 37 R(G) of G in the local Hilbert 
space generated by the Fock configurations T, see Eq. (j3j). 

5 |r} = £i? r , r ( 5 )|r'}. (38) 

r' 
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In the original- original basis [Eq. Q] the invariance con- 
dition 



is equivalent to 



gVg- x =V VjgG 



[X,R(g)}=0 Vj6G. 



(39) 



(40) 



Let us assume that the most general transformation U 
that relates the original and the natural basis 



p 

U\0) = |0> 
c/|r) = |n) 



(41) 



commutes with G, i.e., that 

[C/ ) i?( 5 )] = [C/t,i?( 5 )]=0 VjeG. (42) 

It can be shown, see appendix I A 1[ that Eq. (j32j) is ver- 
ified for every group G whose elements do not mix con- 
figurations that belong to different eigenspaces of the 
number operator N. This assumption is obviously veri- 
fied by every geometry group, but excludes, for instance, 
the particle-hole transformation. A first consequence of 
Eq. (|42|) is that the matrix A has the same form in the 
original-natural and in the original-original representa- 
tion. In fact 

V ee Yl ArHWl = £ (XU) Tn |r>(n| , (43) 

IT' Fn 

and from Eqs. (pf0|) and (fl2|) we have that 

[\U,R(g)]=0 VgeG. (44) 

Notice that from the assumed invariance of |^o) re- 
spect to G we have that, \/g G G, 

Prr> = (*o||r'>(r||*o) 
= (*o||5r'>( 5 r||* > 

= {R\g)PR{g)) TT , . (45) 

Using Eq. (|4"Tj) we can easily express P in terms of P° as 
follows 



p = UP°U* 



(46) 



where 



=<*o||»i>(«i||*o>, (47) 

so that, combining Eq. ([^5]) and Eq. (|^6|) . we obtain that 

P° = (U ] R\g)U) P° (U^R(g)U) . (48) 

Eq. (jlH]) is equivalent, because of Eq. (|42]h to the follow- 
ing invariance relation for P°: 



[P°,R(g)] = VgeG. 



(49) 



From the above considerations and Eq. ([23]) we can 
conclude that <p satisfies the same invariance relation 
of the A coefficients of the Gutzwiller projector in the 
original-original basis, i.e., that 



[<j>,R(g)} = VjeG. 



(50) 



In other words, we have proven that <p has the same form 
of A expressed in the original-original Fock representa- 
tion [Eq. gj)]. 

The set of all the (j) matrices that satisfy Eq. ([50|) is 
a linear space V</>. Consequently, there exists a basis of 
matrices {(/>&} such that 



[<p k ,R(g)}=0 VgeG 

<t> = c fc <p k . 



(51) 
(52) 



In this work we assume that c& and <j)k are real. This 
means that V</> is a linear space over the field of real num- 
bers. Notice that this does not restrict the variational 
freedom as long as the local interaction H\ oc is real. This 
excludes, for example, the spin-orbit coupling. 



A. Calculation of 

In order to calculate {(j)k} it is convenient to apply a 
similarity transformation V to R(G) 



R v (g) = VR(g)V^ Vg € G 



(53) 



with the property to decompose R(G) in irreducible rep- 
resentations #2S More precisely, we need to calculate a uni- 
tary matrix V such that R v (G) is of the form 



R v (g) 



fRY(g) 



R y s {g) J 



V«eG, (54) 



where (i) the representations R% "(G) are irreducible for 
all i G {l,..,s}, (ii) if two representations RY,RY are 
equivalent then they are also equal. In appendix |A] we 
derive a possible procedure to calculate explicitly the sim- 
ilarity transformation V utilized in this section for a gen- 
eral geometry group G. 

Let us consider the linear space Vjf (over the real field) 
of the complex matrices (f) V that satisfy the following 
equation 



[R v {g),4> v }=0 VgeG. 



(55) 



It can be easily proven by means of the Schur lemma^i 
that the most general (j) V G VY is of the form 



(pX 







(56) 
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where the blocks k G {1, .., s} correspond to inequivalent 
representations of G, and each block is of the form 



Pk 



r n k lU k 



r ln k Idfc 



(57) 



td k being identity matrices of size dk x d^, d& being the 
dimension of each one of the irreducible equivalent repre- 
sentations of G repeated in the fc-th block, and r^- being 
independent complex numbers. Eqs. ([56]) and (|57|) allow 
to define straightforwardly a basis {4>X} of V^. 

Let us define now the linear space of real matrices 
generated by the set of matrices {4>k} obtained as 



4> k =vmv. 



(58) 



The linear space that we need, see Eq. (jSUJ), is ob- 
tained as 



v, = v,nw R , 



(59) 



where Wr is the linear space of all real matrices. It is 
convenient, finally, to ortho normalize the basis set 
of Vd> in order to have 



Tr(^) = 6. 



B. Independent variational parameters 



(60) 



All the relevant quantities that define the variational 
energy, i.e., (i) the Gutzwiller constraints [Eqs. ([25i 
126]) ]. (ii) the K matrices [Eq. ([29])] and (hi) the local- 
interaction energy [Eq. (|27j) ] can be expressed in terms of 
quadratic forms in the Ck coefficients defined in Eq. ([52]) 
as follows: 



E 
E 
E 



ap 



1 M^ + M^ 



c), 



(61) 



ij 

= J^CiCjU* = (c\U\c). (63) 

ij 

Notice that the tensors M s ', A/" 5 and J7 are fully de- 
termined by the symmetry of the wavefunction [Eq. (pQ)] 
and the number of orbit als. For this reason it is gener- 
ally convenient to precalculate them before starting the 
numerical minimization of the variational energy. This 
point will be further discussed in Sec. IVII D[ 

C. Simplified variational ansatz 

The complexity of the numerical problem is consid- 
erably reduced if the mixing of different atomic config- 
urations is neglected in the Gutzwiller projector. This 
amounts to assume that the matrix <p defined in Eq. 
has the form 



h 



(64) 
(65) 



where P^ nt are the orthogonal projectors onto the 
eigenspaces of the local atomic interaction Hint- Al- 
though the variational parameters neglected in Eqs. ([6H 
[65]) can play a crucial role in some case^ this simplified 
ansatz merits to be mentioned for at least two reasons, (i) 
It still allows to solve exactly the problem in the atomic 
limit, (ii) The number of independent variational param- 
eters is generally extremely lower in this approximation, 
allowing to perform calculations not feasible otherwise. 
Furthermore, once a variational result 



E c o n £c 1 



(66) 



is obtained assuming Eqs. ([64H65p , it can be used as a 
good starting point c for the self-consistent search of the 
energy minimum with the more general variational space 
discussed before, see Eq. ([52]) , 



c fc = Tr(<^ ) : 



(67) 



with the result to speed up the calculation. Notice, in 
fact, that [H{ nt ,G] = 0, implying that 



[<h,R(g)] = o v^gG, 

i.e., that 0o G V</>. 



(68) 



i° - 



ij 

= E c ^ N % = E \ {*% + O 

ij ij 

= <c|^|c), (62) 



V. NUMERICAL OPTIMIZATION OF THE 
VARIATIONAL ENERGY 

In this section we discuss in detail the self-consistent 
numerical strategy to minimize the energy [Eq. ([27]) ] ful- 
filling the Gutzwiller constraints [Eqs. (j25H26|) ] . 
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Notice that the formulation of the Gutzwiller prob- 
lem through Eqs. ()25ti27|) is formally analog to the con- 
strained formulation of DFT derived by Lev y 38 1 39 and 
Lieb. 40 In fact, the variational energy can be expressed 
as a functional of the variational density matrix n°, see 
Eq. (|36]), 



£var 



[n°] = minfvar [c, * ] , 



(69) 



where min n o denotes the minimum over the set of varia- 
tional parameters c and |^o) 5 satisfying the Gutzwiller 
constraints [Eqs. (|25ti26j) ] at fixed n°. In this work 
Gutzwiller the problem is solved by calculating the den- 
sity functional £ V ar[^°], and minimizing it respect to n°. 

For clarity reasons we have structured the rest of this 
section as an exposition of the numerical procedure, omit- 
ting the mathematical proofs. The mathematical details 
can be found in the appendices. 



A. Preliminary calculation 

Let us consider the Gutzwiller renormalized non-local 
tight binding operator 

T° = E E (70) 

a(3 R/R/ 



where 



^RR' — 2-^ RR' '^ui'^PS • 

a(3 



(71) 



For later convenience, we consider a general one-body 
Hamiltonian 

H G = f G + AH = E 4n vLv kn , (72) 

kn 



where AH is a given local operator. 

Let |^o) be the ground state of H G . It can be easily 
verified that 



B. Slater determinant optimization step 

For later convenience, in this section we solve the prob- 
lem to calculate the state |^o) that realizes the minimum 
of the Gutzwiller variational energy at fixed c 

£ n o, n = min (*|f G |*) 

|*)es„ 

S n o EE {|*> t.c. <*| dl a d R0 |*) = 8 afl n a } , (76) 

where T G is given by Eqs. ([70H7ip . Note that the func- 
tional £ n o^n depends on c only indirectly through 7£, that 
is given by 



(c\M^\c) 



(77) 



see Eqs. ([29]) and (l6Tj) . 

It is convenient to account for the Gutzwiller con- 
straints employing the Lagrange multipliers method. We 
introduce 



H G [R,\] = T G + AH , 



where 



A ^ = Yl X ^ d na d np ■ 

R a(3 



(78) 



(79) 



Notice that H G has the same form of Eq. ([72]) . Finally, 
we calculate the ground state |^o) of H G [JZ, A] for \ a p 
such that the Gutzwiller constraints [Eqs. ((254126)) ] are 
satisfied. Once the Lagrange multipliers \ a p are known 
we compute Eq. ([73]) . 

In summary, the calculations described in this section 
associate the input variables and 1l a p to the output 
matrix 



a/3 



d(y \f G \v ) 

dTlaB 



(80) 



Fig.ffl 



m^MW* (73) 



where 



(/k)nm — @(~ 6< kn) $nm 

VL = E( C/k )«" d ka- 



(74) 
(75) 



Eq ([73]) will be used in the following subsections, where 
two important inner parts of our numerical scheme are 
described in detail. 



C. (j)- matrix optimization step 

In this section we derive the numerical procedure to 
minimize the Gutzwiller energy functional 



£* [c] = (y \f G \y ) + (c\U\c), 



(81) 



where T G is given by Eqs. (|70H7ip and 1Z depends on c 
through Eq. ([77]) . keeping the Slater determinant |^o) 
fixed and respecting the Gutzwiller constraints 



(c|c) = 1 
(c\N^\c) = 5 aP n a . 



(82) 
(83) 
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FIG. 1: (Color online) Flow chart representing the numerical 
calculation of the functional T n o . 



In order to solve this problem we adopt the following 
linearization procedure, that is based on appendix [B] We 
consider the Hermitian matrix 



F[D,X]=H[D]+L[X], 
where V a p is defined by Eq. ([80]) , and 



H[V] = U + Y. V * 



M s 



/ j c*p I 

a/3 y/ n °p( 1 - n p) 



a(3 



(84) 

(85) 
(86) 



Then, we calculate the ground state c of F[D, A] for A 
such that the Gutz wilier constraints [Eqs. ()83ti83|) ] are 
satisfied. The obtained vector c is used to define a new 
1Z through Eq. (I77|) . 

Notice that the matrix V a p entirely encodes the depen- 
dency of the problem on |^o) in Eq. ([84]) . In summary, 
the above calculations associate to the input variables 
and V a p the output renormalizaton matrix 7Z a /3, see 

Fig.m 



We underline that the size of 1Z is equal to the number 
of orbit als, and that the number of independent parame- 
ters that define it is often reduced by symmetry. For this 
reason the solution of Eq. ([55]) generally requires a few 
Newton steps to converge, see Sec. IVII Dl The problem 
of the exponential scaling of the local many-body space 
affects the numerical algorithm exclusively through the 
solution for the ground state of F[D, A], see Eq. (|5^|) . The 
size of the matrix F[D, A] is, in fact, equal to the dimen- 
sion of the vector c. Nevertheless, the calculation of the 
ground state of F[D,\] is not numerically problematic 
for two reasons, (i) The dimension of c is reduced by 
symmetries, as will be shown in Sec. IVII Dl (ii) The cal- 
culation of the ground state of F[D,\] does not require 
a full diagonalization. Less computationally demanding 
algorithms, such as the power method or the Lanczos 
method, can be employed. See Sec. IVII D I for further 
details. 

We remark that the self-consistent numerical algorithm 
derived in this paper requires, as a starting point, only a 
initial "guess" for the variational density matrix n° and 
the matrix 1Z. It is not necessary to construct a good 
initial guess for the whole Gutzwiller wavefunction, i.e., 
for the matrix (/>, while this would be necessary in order to 
perform a direct constrained minimization of the energy 
functional [Eq. ([27]) ] . This implies that the stability of 
the algorithm is not affected by the exponential scaling 
of the number of parameters involved in the calculation. 
This point will be further discussed in Sec. IVII D I 



VI. APPLICATION TO LDA+GUTZWILLER 



D. Fixed point formulation 

A very important observation in our implementation 
is that the composition of the two optimization steps de- 
rived in Sees. IV Bl and IV CI can be described as a func- 
tional T n o that associates a given renormalization matrix 
IZi to a new renormalization matrix IZi+i 

n i+1 =TAKi], (87) 

see Fig.[U This operation lead to a reduction of the varia- 
tional energy unless, by definition, 1Z solves the equation 

T n o[K]-K = 0. (88) 

In this case 1Z defines a stationary point of the energy 
functional. This observation amounts to formulate the 
minimization of the variational energy at fixed n° as a 
fixpoint problem, that can be solved in several ways^i 
A first possibility is to use IZi as an input to obtain 
7Zi+i and iterate the procedure up to convergence. This 
procedure is commonly referred to as forward recursion 
methodA^ 2 - However, as it will be shown in Sec. IVII D] 
the application of the Newton method is generally much 
more efficient. 



In this section we briefly discuss how to combine the 
Gutzwiller scheme with a first principle calculation of the 
uncorrelated electron structure as an input, applying the 
DFT scheme with LDA. This combined scheme is named 
LDA+G. We also outline how the Gutzwiller solver has 
to be modified in order to account for the double count- 
ing^ The double counting appears as a mean-field con- 
tribution in the exchange-correlation taken into account 
in the LDA calculation. 

As a starting point we consider the Kohn-Sham refer- 
ence system obtained within a converged LDA calcula- 
tion 

^DA=EE £ k« S, ?U» ( 89 ) 

k n 

tL|o> = IV>L S >- (90) 

In order to be able to apply LDA+G the tight binding 
Hamiltonian [Eq. ([59]) ] must first be expressed in terms of 
a proper localized basis set 20 \xka)-> where a labels both 
spin a and orbital a. This can be done by means of the 
overlap matrix 

[SkUHV&k.), (9i) 
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giving 



H L BA - 5Z5Z fc k c ka c k/3 
k a(3 

e \f = I^klcm 6 kn [*Sk]n/3 

n 

4a|0) = Wkcc). 



(92) 

(93) 
(94) 



In order to express Hlda in the same form of Eq. (|2]) 
we separate Hlda in a non-local part T and in a local 
part (the crystal fields) as follows 



#lda = T + L 

k a/3 

l = y^ y^ ^ c L c Rg 

R a/3 



where 



.a/3 
k 

_ a/3 _ 7 a/3 
6 k — e k 6 ' 



(95) 
(96) 

(97) 



(98) 
(99) 



and ^ is the number of sites R. The on-site electron 
interaction can be modeled by the Slater-Kanamori ro- 
tationally invariant atomic interaction^ 



(100) 



R 



a^6 era' 



' C Racr C Ra-cr C R6-cr C 



R6cr 



cr 



J 7 



C Rat C Rc4 C R6t C 



R6| ' 



(101) 



We underline that the form [Eq. ()1Q1 j) ] for i7 int is ob- 
tained by implicitly assuming that the single-particle ba- 
sis \xka) is given by real orbitals, i.e., the cubic (crys- 
tal) harmonics. It can be proven that the condition 
U = U' + J + J' ensures the rotational invariance in 
the orbital space. The additional condition J = J' can 
be assumed whenever the spin-orbital coupling is negli- 
gible^ 

Our model is now defined in the form of Eq. (|2]), with 
T given by Eq. (|96|) and H\ oc = Hi nt + L. An additional 
on-site term, the double counting, needs to be added to 
Eq. (|27|) as the average orbital-independent interaction 
energy is already included in LDA. A common choice of 
the double counting term k) 21 i 22 

£dc[p] = ^n(n-l)-£^n CT (rv-l) (102) 



U = 



U + 2U 



21 + 1 

J = tj-U' + J, 



(103) 



where I is the angular momentum quantum number of 
the considered localized basis set, 



n 



(104) 



see Eq. (|3"5|l. and n a is the mean value of Cp^c-p^ with 
respect to the Gutz wilier wavefunction, that is given by 

n a = £c;c/It(</>t/t/ a «/>.) 

= J2 c ^ P «« = (c\ P ™\c). (105) 



The presence of the double counting term gives rise to 
the following additional term in Eq. 



(106) 



When the self-consistent Gutzwiller calculation is con- 
verged the result can be fed back to the LDA code. By 
calculating the density matrix 

Wa^ = (*0| ^4^^1*0), (107) 

and representing it in the Kohn-Sham basis 

[Pk S ]nm = ^2 iSk]na [pk]aj3 [£k]m/3 , (108) 

a(3 

it is possible to get a prescription how to transform the 
Kohn-Sham eigenfunctions and occupancies in order to 
reproduce the physical Gutzwiller electron density^ To 
the new total density corresponds a new effective poten- 
tial for the Kohn-Sham reference system that, in turn, 
defines a new tight binding Hamiltonian. 

The procedure is iterated until self-consistency is 
reached. 

VII. ILLUSTRATIVE RESULTS 

This section has two purposes, (i) As a proof of con- 
cept we present some numerical result in comparison with 
other calculations based on different methods, (ii) We 
outline several technical details of the calculations and 
the speed of convergence of the algorithm. 

Test calculations have been performed on multi-orbital 
model Hamiltonians of the form 



k6cr 



k<7 ab 
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+ EE^R-'W+iJint, (109) 

Rcr ab 

where the structure of H{ nt was defined in Eq. (|101|) . A 
paramagnetic Gutzwiller wavefunction has been assumed 
in all the calculations shown in this section. The hopping 
matrix has been set up as either (i) nearest neighbor hop- 
ping on a three dimensional cubic lattice, giving a non- 
interacting DOS with cusps close to half the bandwidth, 
or (ii) nearest neighbor hopping on a Bethe graph with 
infinite coordination number, corresponding to a semicir- 
cular non- interacting DOS. In both cases the half band- 
width W is set as the unit of energy. 

A. Two-bands Hubbard model 

First of all, let us consider the case of two orbitals. 
In the special case of half-filling and degenerate bands 
we compare our calculations with the available results 
obtained from Ref. |43| by means of the rotationally in- 
variant slave-boson technique^"— that is equivalent to 
the Gutzwiller variational method on the mean field 
level i 48 1 49 In this case ££ 6 is set up as nearest neighbor 
hopping on a three dimensional cubic lattice 

1 3 

with t® b = S ao . For this specific model, in which the 
single-particle energy dispersion of the two bands are 
identical, we have that 

n ab = Vzs ab , (in) 

where Z can be interpreted as a measure of the quasi- 
particle renormalization weight. 50 As shown in Fig. [2j 
the Gutzwiller calculation gives the same values of Z 
as a function of U and J/U as the slave-boson cal- 
culations. The Brinkman-rice transition^ occurs at a 
strongly J-dependent critical U. This well known fact^i 
supports the argument that the spin-exchange on-site in- 
teraction needs to be taken into account to accurately 
describe strongly correlated systems. Furthermore, while 
for J/U = the phase transition is second order, at finite 
J/U it becomes first order, with a hysteresis region char- 
acterized by an additional fixpoint solution of Eq. ([88]) , 
corresponding to a second stationary point of the varia- 
tional energy, see Fig. [2j To our knowledge this solution 
has never been reported before neither in the Gutzwiller 
nor in the slave-boson approximation. 

Let us consider the more general case of non- 
degenerate bands with finite crystal field splitting A 

E /ab 4w<W = A (n Rla - n R2a ) . (112) 

ab 

This model has been studied in detail with DMFT in 
Ref. [52|. Here we compare our Gutzwiller results with 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 

U 



FIG. 2: (Color online) Comparison of results, for two degen- 
erate bands on the 3D cubic lattice with rotational- invariant 
Hunds interaction, from Gutzwiller (solid lines) and slave- 
boson^ (dotted lines), showing the quasi-particle weight Z 
as a function of U and (from right to left), J/U — 0, 0.01, 
0.02, 0.05, 0.10, 0.20, 0.45. Additionally, also the higher en- 
ergy fixpoint solutions of Eq. ((88} at finite J/U are reported. 

part of the available DMFT data. This gives an oppor- 
tunity to discuss some features of the specific implemen- 
tation derived in this work and to introduce some general 
merits and limits of the Gutzwiller variational method in 
itself. 

In Fig. [3] the expectation value of the filling per spin 
hicr is shown for several values of U and J/U. Following 
Sec. IIII| the expectation value of h\ G is given by 

n la = Tr(^flf la cj>), (113) 

see Eq. In the same figure the DMFT results from 

Ref. [52| are also shown. These calculations were per- 
formed assuming a semicircular density of states. 

For the specific crystal field splitting considered, A = 
0.2, the system is metallic at small ?7, and is driven to- 
wards a Mott insulating or an orbitally polarized phase, 
depending on the value of J/U, upon increasing the in- 
teraction strength U. Notice that the DMFT and the 
Gutzwiller results are in very good agreement in the 
metallic phase, especially away from the Mott transition. 
This confirms that the quality of the Gutzwiller calcula- 
tions is generally comparable with the quality of DMFT 
for the ground state properties of strongly correlated met- 
als^ On the contrary, the Mott insulating phase can not 
be described correctly by means of a Gutzwiller approx- 
imation. Nevertheless, it is correct to assume that Z 
approaching indicates that the metallic phase becomes 
unstable, compatibly with the Brinkman-rice scenario^ 
Notice that the critical coupling U C (J) of the Mott tran- 
sition predicted by the Gutzwiller approximation is not 
accurate in general. For instance, at J = the Gutzwiller 
calculation gives U c « 5 for the semicircular density of 
states, which differs from the DMFT result^ [/DMFT _ 4 
by 20% (not shown). Notice that the results of Fig.[2]were 
obtained assuming a three dimensional cubic lattice, and 
not a semicircular density of states. 
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FIG. 3: (Color online) Filling per spin of orbital 1 for 
A = 0.2 and different values of J/U. From bottom to top, 
J/U = 0,0.01,0.02,0.05,0.1,0.15,0.25. The continue lines 
correspond to our Gutzwiller results for the metallic phase 
at half-filling. Open (full) symbols correspond to metallic 
(insulating) solutions obtained from DMFT« at inverse tem- 
perature j3 = 25. 

B. Five-bands Hubbard model 

In this section we study the Hamiltonian of the 
form [Eq. (|109j) ] for five bands, describing correlated d- 
electrons in a cubic crystal. We assume a semicircular 
density of states, and that the full rotational symmetry 
is broken by a finite crystal field splitting A between the 
three ti g and the two e g orbitals, i.e., 

E rbc Ra CT c R6ff =A J2 ( 114 ) 

ab a(zt2 g 

This model has previously been used as a benchmark sys- 
tem in DMFT^ In our calculation we assume a param- 
agnetic Gutzwiller wavefunction invariant with respect to 
the symmetry point group of the cube. This allows to re- 
duce considerably the number of variational parameters, 
see Sec. NTTD\ 

Let us consider the expectation value S 2 of the total 
spin squared S 2 , for which DMFT data are available in 
Ref. [EI From Eq. flU) we have that 

S 2 = Tr(tf S 2 (j>) , (115) 

where 

and a k are the Pauli matrices. In Fig. H] the behavior 
of S 2 is shown at fixed U = 1 for 6 electrons per site 
and several values of J/U. In the figure our results are 
compared with the DMFT data from Ref. [H. Consis- 
tently with DMFT, we find that S 2 grows monotonically 



FIG. 4: (Color online) Gutzwiller expectation value of S 2 for 
different values of J/U and crystal field splittings A, the total 
filling is 6 electrons per site and U = 1. Comparison with 
DMFT data at inverse temperature f3 = 25 from Ref. [53|. 

upon increasing J/U, and that the crystal field splitting 
A = 0.25 slightly reduces S 2 compared to the case of de- 
generate bands. Notice that the discrepancy between the 
Gutzwiller results and the DMFT data becomes larger 
upon increasing J/U at fixed U. A similar qualitative 
behavior could be observed even in the calculation shown 
in Fig. [2] for the two-bands model. Nevertheless, the ob- 
served deviation between the Gutzwiller results and the 
DMFT data seems to be more substantial in this case. 



C. Bilayer Hubbard model 

In both the models previously considered the renor- 
malization matrix 7Z was diagonal due to symmetry. 
For completeness, we also consider the bilayer Hubbard 
model j 32 i 34 i 43 in which 7Z have finite off-diagonal ele- 
ments. In particular, we consider the Hamiltonian given 
by Eq. (|109p assuming that the local hybridization term 
described by the matrix 

>={°vl) <-> 

with V = 0.25, and a hopping matrix given by 
Eq. (|110p as in Sec. IVII A[ Finally, we assume that 
the local interaction Hi nt is given by Eq. (|10ip with 
U' = J' = J = 0. This model has previously been 
studied - with the same parameters - in Ref. 43 with 
the slave-boson method. 

When the bandwidths are equal for the two bands (as 
in the present case) the matrix I defined in Eq. (|117p can 
be diagonalized without modifying the hopping matrix 
for any ki^ 3 - This change of basis transforms the hy- 
bridization term I in a crystal field splitting between the 
bonding (+) and antibonding (-) orbitals. In the new ba- 
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FIG. 5: (Color online) Gutzwiller renormalization matrix Z 
and filling of the bonding-antibonding bands for the two- 
bands bilayer Hubbard model with equal bandwidths, local 
hybridization V = 0.25, U' = J' = J = 0, and filling 
N = 1.88 (solid lines). Comparison with slave-boson results 
from Ref . \M (dotted lines) . 



sis both the effective renormalization matrix Z° = [JZ°] 2 
and the density matrix are diagonal 



Z- 



n+ 
n_ 



(118) 



and the coefficients Z+ , Z- can be interpreted as the 
quasi-particle renormalization weights of the bonding 
and antibonding orbitals respectively. In the original ba- 
sis, instead, 7Z has non-zero off-diagonal elements and 
Z = 1Z 2 has the form Z\\ = Z22 and Z\2 = ^21, where 



2n 
I^i2| 



2 

|Z+ - 2_ 



(119) 
(120) 
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FIG. 6: (Color online) Convergence of the forward recursion 
scheme (circles) and of the Newton method (squares) for 5 
bands at TV = 5, J = and A = 0. The quasi-particle weight 
Z = 7Za a goes to zero as U goes to the critical coupling U c « 
10 . Simultaneously, the leading eigenvalue A of the Jacobian 
of the recursion function goes to 1, and the number of forward 
iterations required to reach a fixed relative precision (here 
\lZi — Hi+i\ < 10 -6 is used) diverges. On the contrary, the 
number of Newton iterations is almost independent of U. The 
matrix [7^o]a/3 = is used as initial condition of both the 
forward recursion series and the Newton method V77. 



In order to compare with the slave-boson results of 
Ref. E3I we have studied the system for TV = 1.88 electrons 
per site. As seen in Fig. [5j the average of Z + and 
given by 2n, decreases monotonically as a function of 
U as expected. Concomitantly, the difference between 
Z+ and Z-^ given by \Z\^\, progressively increases with 
U. Our calculations compare well with the slave-boson 
results, although we find a slightly lower renormalization 
of the antibonding state at large interactions. 



D. Technical remarks 

In this section we point out several technical details of 
the numerical simulations performed to derive the results 
presented above. 

The main technical problem of the Gutzwiller method 
is that the dimension of c, see Eq. (|52]h scales expo- 
nentially with the number of correlated orbitals. Fortu- 
nately, this number can be highly reduced by taking into 
account the symmetries of the system. As an example, 



the number of matrix elements of (j) (i.e., the dimension 
squared of the local Fock space) is compared in table U 
with the dimension of the vector c in the case of a para- 
magnetic Gutzwiller wavefunction invariant with respect 
to the point symmetry group of the cube. This simpli- 
fication is very important, as dim(c) is equal to the size 
of F[D, A], see Eq. J53J), whose ground state needs to be 
evaluated many times during the calculations. To com- 
pute the ground state of F[D, A], i.e., its eigenvector with 
the lowest eigenvalue, we have used the iterative Arnoldi 
based solver provided by the ARPACK library. This cal- 
culation is further speeded up by exploiting the sparsity 
of F[D, A], effectively reducing the cost of the necessary 
matrix- vector multiplications. 



As anticipated in Sec. IIVB( it is generally convenient 
to precalculate <pk and the tensors M^, N^a and U 1 ^ in 
order to further speed up the calculations. This reduces 
the construction of the matrix F to the sums in Eqs. ([85l 
, but increases the memory requirements. In fact, the 
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total number of elements Nt in the tensors scales as, 

N T = (2iV 2 rb + l)iV c 2 + N C N* , (121) 

where N c is the dimension of the vector c, N or \) is the 
number of orbitals, and Nr = 2 2Norb is the dimension 
of the local space. However, the number of stored ma- 
trix elements can considerably reduced by exploiting the 
sparsity of the tensors. For instance, the number of non- 
zero elements in the tensors was reduced by around three 
orders of magnitude for the five-bands Hubbard model of 
Sec. I VII B I Eventually, in more complicated calculations 
it may happen that the number of variational parameters 
is so large that the tensors can not be stored in memory. 
In this case it is still possible to calculate the traces "on 
the fly" . This operation is trivially parallelizable. 

It is well known that the speed of convergence of the 
forward recursion method [Eq. ([87]) ] is limited by the 
magnitude of the largest eigenvalue A of the Jacobian 
of the transformation T n o in the fixed point 1Z. In partic- 
ular, if | A | — ^ 1 the rate of convergence displays a critical 
slowing down. We have found that this situation ac- 
tually occurs in our simulation when U approaches the 
Brinkman-Rice critical value. This is shown in Fig. [6j 
where the convergence of 1Z is shown for different values 
of U at half-filling (N = 5), J = and A = 0. The value 
of A was obtained as 

A = lim (lZ i+1 - K)/(Ki - K) 
Tl = lim Ki , (122) 

where IZi was obtained from the forward recursion se- 
ries [Eq. ([87)) ] starting from the initial condition 

[Ro\ a p = S a p W. (123) 

In DMFT the self energy E is obtained as the solution 
of a fixed point problem^ analogously to the matrix 7£, 
that is the solution of Eq. ([88]) . It is known that also in 
DMFT the rate of convergence of the forward recursion 
method displays a critical slowing down in the vicinity 
of the Mott transition. This behavior has recently been 
shown to be related to the fact that the maximum eigen- 
value of the Jacobian A approaches 1 as U approaches 
the critical valued, as in our Gutz wilier calculations. In 
this case the convergence problem has been cured by em- 
ploying quasi-Newton methods instead of the forward re- 
cursion scheme j 54 i 55 The same strategy is applicable to 
solve Eq. ([88]) . As shown in Fig. [6j this strategy is very 
efficient. While the forward recursion algorithm slows 
down as U approaches its critical value, the number of 
required Newton iterations is essentially independent of 
U. The time required to calculate the results shown in 
Fig. [6] with the Newton method is less than one minute 
for every single U. Nevertheless, the forward recursion 
method could be more stable in some case, as every for- 
ward recursion step leads to a decrease in total energy. 
In fact, see Sec. El every evaluation of T n o corresponds to 
a minimization of the energy with respect to the Slater 



TABLE I: Dimensions of the variational space for 1, 3 and 5 
atomic orbitals (corresponding to s, p and d electrons). The 
number of matrix elements of 0, 2 4iY ° rb , is compared with the 
dimension of the vector c, that is reduced by the point group 
symmetry of the cubic lattice. 





I 


iVorb 


2 4iV orb 


dim(c) 


s 





1 


16 


3 


p 


1 


3 


4096 


16 


d 


2 


5 


1048576 


873 



determinant followed by a minimization with respect to 
the Gutzwiller projector. This guarantees that the fixed 
points calculated by the forward-recursion method are 
local minima, while the Newton method can converge 
also to fix points with one or more Jacobian-eigenvalues 
|A| > 1, i.e., to stationary points of the energy that are 
not local minima. 

We remark that the numerical procedure proposed in 
this paper is divided in two steps, (i) Construction of 
the functional E var [n°] optimizing the variational energy 
for a fixed variational density matrix n°. This optimiza- 
tion can be reduced to the fixed point problem [Eq. ([88]) ] 
and solved with the methods discussed above, (ii) Direct 
minimization of £ V ar[^°] with respect to n°. For com- 
pleteness, this procedure is illustrated explicitly here for 
the bilayer Hubbard model discussed in Sec. I VII Ci For 
each value of U the total energy £ V ar[^°] was optimized 
with respect to the variational density n° using a bound 
minimization routine. An example is shown in Fig. [71 
where the renormalization factors and the total energy 
are shown for fixed U = 2.5 and total filling per site 
N = 1.88 as a function of the antibonding orbital fill- 
ing n_. For each n_ the fixpoint problem of Eq. ([88]) 
was solved employing a quasi-Newton method with con- 
vergence criterion \lZi — 72^+ 1| < 10 -12 in less than 40 
steps. 

Finally, we point out that the presence of finite off- 
diagonal terms in 1Z is due to the generality of the vari- 
ational ansatz considered in this work. In fact, 1Z is al- 
ways diagonal if Eqs. (j30H3ip hold, as it was assumed in 
Ref. [H see Sec. HITCl 



VIII. CONCLUSIONS 

In this article we have derived a numerically efficient 
self-consistent implementation of the Gutzwiller varia- 
tional method. The method proposed was obtained as 
a combination of the self-consistent numerical procedure 
recently derived by Deng et at. in Ref. [20| and the mathe- 
matical formulation of the Gutzwiller problem developed 
by Fabrizio and collaborators ! 30 ! 32 " — This formalism al- 
lows us to overcome the restriction to density-density 
interaction, that was assumed in Ref. [20|, without in- 
creasing the complexity of the numerical algorithm. The 
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FIG. 7: (Color online) Sweep in the filling of the antibond- 
ing orbital n_ for the two-bands bilayer Hubbard model 
with equal bandwidths, local hybridization V = 0.25, filling 
N = 1.88, U' = J' = J = and U = 2.5. Renormalization 
matrix Z (top panel) and total energy (bottom panel). 



approach drastically reduces the problem of the high- 
dimensional Gutzwiller minimization by mapping it to 
a minimization only in the variational density matrix, 
in the spirit of the Lev y 38 1 39 and Liet4£ formulation of 
DFT. For fixed density the Gutzwiller renormalization 
matrix is determined as a fixpoint of a proper functional 
of 7£, whose evaluation only requires ground-state calcu- 
lations of matrices defined in the Gutzwiller variational 
space. We have compared different methods to solve the 
fixpoint problem, finding that the Newton method is gen- 
erally more efficient than the forward iteration method. 
The formalism also allows us to reduce the number of 
independent variational parameters in a well controlled 
way using symmetries. As a proof of concept we have 
performed a few numerical calculations for two and five 
band Hubbard models with full rotationally invariant in- 
teraction, finding good agreement with available DMFT 
and slave-boson data. This analysis shows that the nu- 
merical approach derived is very stable and efficient. For 
these reasons the scheme is promising for first-principle 
studies of real materials, e.g., in combination with DFT 
(LDA+G). 

It is noteworthy that the numerical implementation 
presented in this work allows for straightforward exten- 
sions in two directions of interest, (i) The variational 
freedom of the Gutzwiller wavefunction can be general- 
ized in order to describe superconductin g 32 ! 34 and mag- 
netic systems^ 3 - (ii) The assumption that the coefficients 
A of the Gutzwiller projector [Eq. (0J)] are real can be 
dropped, and generalized to complex values. This allows, 
for instance, to account also for spin-orbit corrections to 
the on-site interaction. 
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Appendix A: Irreducible representation for a 
paramagnetic wavefunction 

In this appendix we explain how to calculate the trans- 
formation V introduced in Eq. ([53]) for a general group 
G. Let us consider the example in which G is the symme- 
try group of a paramagnetic wavefunction invariant with 
respect to a specific discrete group of real (orbital) rota- 
tions G rb, e -&-> the symmetry group of the cube. The 
problem consists in the definition of the most general (j) 
matrix that commutes with the number operator TV, the 
representation of the spin S and the representation G or b 

of Corb- 
ie, TV] = [</>, S] = [0, g] = Vge G orh . (Al) 

In order to achieve our purpose the first step is to di- 
agonalize simultaneously TV and S 2 . Each simultaneous 
eigenspace of these operators is the basis of a represen- 
tation of the symmetry group G identified by the eigen- 
values (TV, £). We denote such a space Vn,s- Let us de- 
compose Vn,s m irreducible representations of the spin 
rotations. In order to do this we consider the kernel of 
the spin lowering operator S- in Vat, 5, and denote it 
by Vn,s,-s- Then, we calculate an orthonormal basis of 
Vn,s,-s m which, for later convenience, also L 2 and Lz 
are diagonal 

V N ,s,s = Span ({^',5,-s}) ■ ( A2 ) 

To each value of L, ttil and i corresponds, by applying 
to s -s the raising S+ operator up to 2s times, a set 

of states labeled as {^'™ L m l g }. Each subset V^™ L ' 1 is 
defined as 

= Span ({Vfei I ™* = - 5 < - S }) ' ( A3 ) 

is a basis of an irreducible representation of the spin 
group. 

It is clear that, in order to commute with TV and S 2 , 
(j) is decomposed in uncoupled blocks, each of them act- 
ing on the corresponding subspace Vn, s- In each block 
we group together the vectors {V^'^ras} w ith the same 
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ms- The Schur lemma 37 ensures that, if the above order 
convention is used, the (TV, S) block of <p has the general 
form 



V 



N,S 







01 



v N „ 



(A4) 



••• p 



N,S 



The structure of the matrix p in Eq. (|A4|) is further 
reduced by the condition [<p : G or b] = 0. The vector space 
V N ,s,m s generated by {^N^s^m 3 }L,m L ,i is the basis of a 
representation of G rb, that can be decomposed in irre- 
ducible representations using standard methods. Notice, 
in fact, that a state transforms exactly as the 

spherical harmonic function under rotations. Each 
one of the obtained irreducible representation of G or b is 
labeled by its characters^ 7 - We group together all the 
equivalent representations with equal characters x- This 
amounts to express VN,s,m s as the direct sum of V$ s ms . 
The Schur lemma 37 ensures that each one of the p N,s 
blocks defined in Eq. (|A4|) has the general form 



\ 







An rh 



(A5) 



where n c h is the number of inequivalent representation 

in VN,S,m S ' 

The final step is to identify the n& states of each sub- 
space Vjy 5 ms that belong to the same row^ of the corre- 
sponding (equivalent) irreducible representations of G or b- 
The Schur lemma 3 - 7 - restricts the structure of each q^ s 
block of Eq. rtA5]) as in Eq. (f57|). i.e., 



( 



N,S, X k- 
11 



N,S, X k 



'ln k ^d k 



r N,S, X kl, I 
num. -"-dfc / 



(A6) 



where td k are identity matrices of size dk x dfc and is 
the dimension of each one of the irreducible equivalent 
representations of G or b repeated in V^ k s . 



ent eigenspaces of the number operator N . We need two 
preliminary observations, (i) By assumption our uncorre- 
cted wavefunction |^o) is invariant respect to the action 
of G, i.e., 



5|*o)=e i *»|* > V 9 GG. 



(A7) 



(ii) From Eq. (j37]) we have that 



where 



D\g)p°D(g)=p° VjeG, 



ftp = <*o| |tf > 



(A8) 



(A9) 



is the variational density matrix expressed in the original 
basis, see Eq. JTDJ). The density matrix p° has exactly the 
same form of Eqs. (|A4IIA6p because of Eq. (|A8|) . For this 
reason it can be diagonalized by means of a matrix U of 
the same form, i.e., 



D\g)UD{g)=U V 5 eG. 



(A10) 



The single particle transformations induced by the matri- 
ces U and D(g) into the Fock many body space evidently 
commute as a consequence of Eq. (jAlQj) . This concludes 
the proof of Eq. 



Appendix B: Self-consistent formulation of the 
0-matrix optimization 



We need to minimize the energy functional defined in 
Eq. (|5Ij) respect to the vector c fulfilling the Gutzwiller 
constraints [Eqs. (|82ti83|) ] . The Gutzwiller constraints 
can be ensured by means of the following Lagrange func- 
tional 



1. Proof of the assumption [Eq. (g2J] 



C[c,X] 



:<w<c|jv^ic). 



a/3 



(Bl) 



Let us prove that Eq. (|^2j) is verified for every group 
G that does not mix configurations belonging to differ- 



The variation of £^ [c] is given by 



S£[c] 



WE 

a(3 



d{^ \f G [<i ,c}\^ } 



+ (5c\U \c) + {c\U \Sc) , 



-n°) 



|C) + (C|E ^^M 



y/n$(l-n$) 



\Sc) 



(B2) 
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and the variation of the Lagrange functional [Eq. (|B1|) ] is equivalent to the following "nonlinear eigenvalue prob- 
is given by lem" 

5C[c,X] = {5c\Y,KpN s ae \c) 

a/3 

+ (c\ KpN^ \Sc) . (B3) F* [c, A] |c> = E\c) , (B5) 

a(3 

The condition that the variation of 
&F* [c, A] = SSy [c] + SC[c, A] = V£c _L c (B4) where 
I 



it r- AW/ I V ^o|r G [^o,c]|^ ) ^ 

A J - ^ + 2^ ™? / = + 2^ • 

y/n°(l-n°) tfi 



(B6) 



In principle, the minimization of £^ [c] could be per- 
formed recursively, starting from a given "guess" Co and 
iterating the following eigenvalue problem 



[ c n, A n ] |c n+ i) — £^ n+ i|c n+ i), 



(B7) 



where £? n +i is the lowest eigenvalue of i^ [c n ,A n ] and 
A n are the Lagrange multipliers such that c n+ i satisfies 
the Gutzwiller constraints. The minimum of £# [c] is 
realized in 



lim c n 

n— ?-oo 



(B8) 



Instead to calculate c m i n , it is convenient to approxi- 
mate c m i n with c\ . This approximation can be considered 
as the result of a "linearization" of the functional £^ [c] 
around the initial guess cq. 



Appendix C: Numerical implementation of the tight 
binding problem 

Let us consider a general translational invariant non- 
interacting tight binding Hamiltonian 



where 



H = T + AH (CI) 



T = E E *RR' (02) 

a(3 R/R' 



ah = E^E^L^- 

a(3 R 

The translational invariance of the system 

,otj3 ,a,f3 W "D R 

RR' ~ r R+R R / +R V±1 0, OL, p 

allows to express H in /c-space 



(C3) 



k/3 ' 



(C4) 
(C5) 



a/3 k 



where 



(C6) 



R 



From Eq. (jC5|) can be easily diagonalized numerically 
and expressed in terms of it's eigenoperators as follows 



(C7) 



kn 

Notice that the overlap matrix 

(C/kU = (0\d ka vL 10) 
appears in Eq. ([75]) . 



(C8) 
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